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ABSTRACT 

Background: The purpose of any numerical technique is to give a result that moves progressively closer to the 
exact solution after each successive iteration. 

Methods: In this paper, we give the exact solution and the numerical approximation of the problem using Euler 
method coded with Maple software to show that the divergence reported earlier can be overcome. 

Results: Presently, the numerical method is convergent for the system of ordinary differential equation using 
Euler method with algorithm written with Maple software even at the region far beyond earlier reported. 

Conclusion: The problem of divergence reported earlier is overcome by using a more appropriate software and 
approach. Convergence and consistence of a numerical solution also depend on the software used to solve the problem. 
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INTRODUCTION 
Background 

Differential equations are used to model problems involving the change of some variable with respect to another. 
Ordinary differential equations are equations in which there is a single independent variable and one or more dependent 
variable. An initial value problem is a differential equation that satisfies a given initial condition. These kind of problems 
occur frequently in application. Their numerical solution is a central task in all simulation environment. 

Numerical methods have been around for a long time. The popular Euler method published in 1768 is attributed to 
Leonhard Euler (1707 - 1783). The method is one of the simplest and most elementary method of solving single and 
systems of ordinary differential equations. 

The application of numerical analysis and numerical methods occurs in just about every field of science and 
engineering. The usage of numerical methods were limited before due to lengthy hand calculations involved in their 
implementation. 

Nowadays, this is no longer the case due to the rapidly changing digital computer industry. Digital computers 
have provided a fast computational device for the development and implementation of numerical methods which can 
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handle a variety of difficult mathematical problems. The advent of different computer software and programming 
languages have made solving different types of various problems easier. 

A numerical method is said to be convergent if the numerical solution approaches the exact solution as the step 
size goes to zero, otherwise it diverges. 

A method might converge on a problem but diverge on another, or converge with one set of starting values but not 
on another. It is essential that a method converges on all problems in a reasonably large class with all reasonable starting 
values for consistency. The divergence and general nature of the numerical solution in Table 2 is the motivation for the 
study. 

It was reported in [13] that increasing the number of steps with Matlab makes the Euler solution drift away from 
the solution. This was what happened with the Euler solution reported in [2]. We seek the solution to the problem with 
another software package which is Maple. 

2.0 METHODS 

In this section, the description of the methods used to arrive at the exact and numerical solutions were explained. 

2.1 The Exact Solution 

Exact solutions are closed-form or implicit analytical expressions that satisfy the given problem. 

In practice, few of the problems originating from the study of physical phenomenon can be solved exactly. 

The system of ordinary differential equations (1) and (2) given in [2] can be solved analytically to obtain the exact 
solution since the equations are linear. Hence we proceed to give the exact solution. 



The model is as follows: 

C'(t) = -(k + 0) C(t)+X N(t) C(0) = Co ( 1 ) 

N'(t) = I(t) - (v + 5) C(t) - |iN(t) N(0)=N 0 (2) 
Where 

9 = y + |i+v + 5 (3) 

N(t) = D(t) + C(t) (4) 

N'(t) = D'(t) + C'(t) (5) 

D'(t) = I(t)-(X+p)D(t) + yC(t) (6) 

C'(t) = I(t) + W(t)- (v + 5 + p + y) C(t) (7) 
C(t) =Number of diabetics with complications 



N(t) =Number of diabetics with and without complications 
I = I(t) = Incidence rate of diabetes. 
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Table 1: Definition of Parameters and Hypothetical Values 



Parameter 


Definition of Parameter 


Value 


s 


mortality rate due to complications 


0.05 


X 


probability of developing a complication 


0.66 


V 


rate at which patients with complications 
become severely disabled 


0.05 


r 


rate at which complications are cured 


0.08 


n 


natural mortality rate 


0.02 


I 


Incidence of diabetes mellitus 


6x 10 6 



Differentiating equation (1) gives 
C"(t) = - ( X + 9) C'(t) + X N'(t) 

Substituting equations (1), (2) into equation (8) result to 
C"(t) = [ ( X + 9) 2 - Uy + 5)] C(t) -XCk + Q + n) N(t) + X I(t) 
From equation ( 1 ) 

N(t) = i(C'(t)-I(t) + a+9) C(t)) 

Substituting equation (10) into equation (9) yields 

C"(t) = - ( X + 9 + u)C'(t) + X I(t) + [( X + 9) 2 - X(v + 5) - ( X + 9 + u) ( X + 9) ] C(t) 

C"(t) = - ( X + 9 + u)C'(t) + X I(t) - ( Xv + X8 + X\x + ^9) C(t) 

Hence we obtained 

C"(t) + oC'(t) + pC(t) = X I(t) 

where 

o = X+ 9 + |i 

P = Xv + Xd + X\i + (j.9 

The auxiliary equation for the homogenous part of equation (13) is 
m 2 + om + P = 0 

Solving the quadratic equation (1 1), we have 

-a±Vo 2 -4R 

m = 

2 

The complimentary solution of equation (8) is 

C(t) = K 1 e" 1/2(a "V oA2 - 4 P) t + K 2 e" 1/2 (o + V 0 " 2 - 4 ?) 1 

Using the method of undetermined coefficient to find the particular solution 
Let 
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Cp(t) = K (19) 

C'p(t) = 0 (20) 

C" P (t) = 0 (21) 

Substituting equations (19), (20) and (21) into equation (13) we obtained 

Cp(t) = ^I(t) (22) 
The complete solution of equation (13) is 

C(t) = Kje-^^-V^^Dt + K2 e-i/2 (a + 1 /^4p)t + jj I(t) (23) 
Let 

ti! = i(o-Vo A 2 - 4P) (24) 

tl2 = ^(o + V oA2 - 4 P) (25) 
C(t) = K^-i 1 ' + K 2 e-i 2 ' + jj I(t) (26) 

C'(t) = - ( iiiKie" 1 ! 1 1 + r) 2 K 2 e-V ') (27) 
Substituting equations (21) and (22) into equation (5), we have 

N(t) = K^-i 1 1 + K 2 e-i 2 '+ - I(t) + -Kie-i 1 1 + -K 2 e^ 2 1 + - I(t) 

ft X x p 

-idl.K^-i 1 ' +i! 2 K 2 e-i 2t ) (28) 

Using the initial conditions on equations (21) and (23) and simplifying, we obtained 

p(X+e-ii2)C0+ XI(t)ii2 -XpNO , om 

P0ii-n2) 

-p(X + 9-til)C0-U(t)til+XpN0 ,„ m 

JS. 2 = (jU) 

P (11-12) 

3. 0 RESULTS AND DISCUSSIONS 

In this section we discuss the results obtained using Mable software. 
3. 1 Exact Solution 

Substituting the values of each parameter as given in Table 1. The exact solution of the problem is as follows: 
-^j (-11+/69")'/ 19300OOO r77T 3875000 



-— (u+JW)t, 1930 0000 r— 3875000 
+ £ [ 299 V69-— — 

618750000 
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N ^ 7 e lJ 19300000 jog 3875000' 



11 I, 299 13 ) 

2 (-11+769")*/ 19300000 ,-t- 

+ 3T { 2^^ 69 

3875000 ^ ,— , 7 -jj-Ch +■/»")'/ 19300000 r— - 
1 ' 69 + — — e — — J 69 



■)- 



13 J v 11 I 299 

3875000 ^ 2 "25" (" +^"> ' { 19300000 



" e " " " V 69 



13 J 33 I 299 



3875000 ^ r-rzr , 806250000 
— 13— JV» + [3 



3. 2 The Numerical Solution 

In this section, we discuss the numerical methods used and the results obtained with different software packages. 
3. 2. 1 The Euler Method 

The method is applied as follows: 

C n+1 = C n + h(-(X + 9)C n + XN n ) (31) 
N n+1 = N„ + h ( I- (v + 5) C„ - uN n ) (32) 
where h is the step size. The sequence is constructed as t 0 , ti = t 0 + h, t 2 = t 0 + 2h, . . . 

The numerical solution for the system of equation (1) and (2) using Euler method and Pade approximation 
(Method 1) coded with Matlab as reported in [2] is given below for the purpose of comparison and illustration. 



Table 2: Table of Numerical Solution with Matlab 





N(t) 


N(t) 


C(t) 


C(t) 


t 


Euler 


Method 1 


Euler 


Method 1 


0.01 


6.08 


6.08 


4.75 


4.75 


0.02 


6.11 


6.11 


4.77 


4.77 


0.05 


6.11 


6.11 


4.7 


4.7 


0.1 


6.11 


6.11 


4.7 


4.7 


0.2 


6.11 


6.11 


4.7 


4.7 


0.5 


6.11 


6.11 


4.7 


4.7 


1 


6.11 


6.11 


4.7 


4.7 


2 


6.11 


6.11 


4.7 


4.7 


2.5 


div 


6.11 


div 


4.7 


3 


div 


6.11 


div 


4.7 


3.5 


div 


6.11 


div 


4.7 


4 


div 


6.11 


div 


4.7 



div = divergence of numerical method. 



Table 2 shows that the result for Euler solution converges until 2.0 and diverges thereafter. The solution were also 
given to 2 decimal places for N(t) and 1 decimal place for most of C(t). This doesn't depict a good approximation. 
The solutions as discussed in [2] for Euler and Pade approximation from the second iteration to the end were all constant. 
This also is not the best for an approximate solution. 
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The numerical solution for system (1) and (2) using Euler method is given below with the exact solution. 
The algorithms were coded in Maple 14 software. 



Table 3: Numerical and Analytical Solution with Maple Software 





N(t) 


N(t) 


C(t) 


C(t) 












0.01 


6.110078000 


6.110078039 


4.699906000 


4.699906660 


0.02 


6.110156078 


6.110156155 


4.699813323 


4.699814632 


0.05 


6.110390770 


6.110390952 


4.699543124 


4.699546316 


0.1 


6.110783369 


6.110783702 


4.699118258 


4.699124385 


0.2 


6.111573422 


6.111573967 


4.698359425 


4.698370704 


0.5 


6.113972187 


6.113972793 


4.696726193 


4.696748133 


1 


6.118009418 


6.118008833 


4.695681305 


4.695709876 


2 


6.125946876 


6.125942007 


4.697419581 


4.697442136 


2.5 


6.129752812 


6.129745688 


4.699324926 


4.699341454 


3 


6.133417059 


6.133407856 


4.701566698 


4.701577185 


3.5 


6.136928412 


6.136917366 


4.703987005 


4.703992022 


4 


6.140282245 


6.140269605 


4.706481314 


4.706481643 


4.5 


6.143478236 


6.143464243 


4.708981156 


4.708977607 


5 


6.146518841 


6.146503716 


4.711442411 


4.711435728 


5.5 


6.149408258 


6.149392200 


4.713837350 


4.713828185 


6 


6.152151739 


6.152134919 


4.716149268 


4.716138162 


6.5 


6.154755112 


6.154737687 


4.718368818 


4.718356215 


7 


6.157224497 


6.157206599 


4.720491546 


4.720477807 


7.5 


6.159566085 


6.159547830 


4.722516221 


4.722501635 


8 


6.161786011 


6.161767501 


4.724443707 


4.724428502 


! 8.5 


6.163890279 


6.163871596 


4.726276197 


4.726260555 


9 


6.165884685 


6.165865910 


4.728016703 


4.728000769 


9.5 


6.167774822 


6.167756021 


4.729668713 


4.729652602 


10 


6.169566041 


6.169547270 


4.731235956 


4.731219768 



Table 3 shows the results of the exact solution and numerical solution respectively using Euler method coded with 
Maple 14 software. The results were given to ten significant figures which is appropriate for an approximation method. 
The result converges at all point from 0 to 10, a region beyond the one in [2]. The results also depict iterative nature of 
numerical methods. 

The error which is the difference between the exact solution and the numerical approximation at specified time 
values is relatively small. The Euler method gives a reasonably good approximation to the exact solution. 

4. 0 CONCLUSIONS 

In this present research study, we have shown that the problem converges continuously even at a far higher region 
with Euler method coded using Maple. The computed values are close to the analytic solution. The method is constructive, 
effective and reliable using Maple compared to Matlab. Numerical instability and divergence can also depend on the 
software used. A necessary and better approach could be to always compare with the exact solution where applicable 
especially for linear problems like the present problem and to verify with other software or programming languages. 
However no code is infallible, and when a spurious result is obtained from a code the data can be rationalize with the aid of 
the code's documentation as backup for the numerical solution. 
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